Discrete soliton ratchets driven by biharmonic fields 
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Directed motion of topological solitons (kinks or antikinks) in the damped and AC-driven discrete 
sine-Gordon system is investigated. We show that if the driving field breaks certain time-space 
symmetries, the soliton can perform unidirectional motion. The phenomenon resembles the well 
known effects of ratchet transport and nonlinear harmonic mixing. Direction of the motion and its 
velocity depends on the shape of the AC drive. Necessary conditions for the occurrence of the effect 
are formulated. In comparison with the previously studied continuum case, the discrete case shows 
a number of new features: non-zero depinning threshold for the driving amplitude, locking to the 
rational fractions of the driving frequency, and diffusive ratchet motion in the case of weak intersite 
coupling. 
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I. INTRODUCTION 

Transport phenomena induced by the interplay be- 
tween non-equilibrium fluctuations, symmetry breaking 
and nonlinearity, have recently attracted a great deal of 
interest. In particular, point particle ratchets described 
by ordinary differential equations have been largely inves- 
tigated due to their relevance in several fields, including 
molecular motors and Joscphson junctions (see reviews 
[1,2]). In simple terms, the point particle ratchets ap- 
pear as the unidirectional motion of a damped and driven 
particle, which is achieved under the influence of only 
stochastic and deterministic forces of zero average, in- 
dependently on initial conditions. The phenomenon was 
ascribed to the breaking of the symmetries connecting 
orbits with opposite velocities in the phase space [3-5] 
and to the phase locking of the particle dynamics to the 
external driver [6,7]. 

Ratchet phenomena in infinite-dimensional systems de- 
scribed by nonlinear partial differential equations of soli- 
ton type have also been investigated. This has been done 
for both spatially asymmetric potentials with temporar- 
ily symmetric AC fields [8,9] and for symmetric poten- 
tials with temporarily asymmetric AC fields [10-13]. In 
both cases, the ratchet effect appears as a unidirectional 
motion of the soliton, which resembles the drift dynam- 
ics observed for point particle ratchets; from here the 
name of soliton ratchets originates [9] . The soliton ratch- 
ets induced by asymmetrical external fields have been 
implemented experimentally in long Josephson junctions 
[14] by means of nonlinear harmonic mixings. This ap- 
proach has been shown to be effective also in other phys- 
ical contexts [15] (for a detailed list of references see 
the review paper [2]). The existence of soliton ratch- 
ets in long Josephson junctions was also experimentally 



demonstrated by using asymmetric magnetic fields [16] 
and spatially asymmetric currents [17]. From the theo- 
retical point of view, symmetry breaking conditions to 
generate soliton ratchets were discussed in Ref. [11]. The 
mechanism underlying soliton ratchets was proposed in 
Ref. [9] for the case of a perturbed asymmetric double 
sine-Gordon equation driven by a symmetric AC driver 
and extended in Ref. [10] to the case of the damped sine- 
Gordon (SG) equation with asymmetric AC fields. In 
both cases, the phenomenon was ascribed to the exis- 
tence of an internal oscillation on the kink profile which, 
in the presence of damping, couples to the translational 
mode of the kink and produces transport. The internal 
vibration was shown to be spatially asymmetric, thus giv- 
ing directionality to the motion, and phase locked to the 
external force. This mechanism, also known as the inter- 
nal mode mechanism, has been confirmed for a number 
of systems such as the asymmetric double sine- Gordon 
equation with symmetric driver [9,18] and the SG sys- 
tem with temporal asymmetric forces [10,19,20]. 

In contrast to the continuous case, however, discrete 
soliton ratchets have been scarcely investigated (some 
work on spatially asymmetric discrete soliton ratchets 
has been done in Refs. [21,22]). In this case, one can ex- 
pect that the presence of the Peierls-Nabarro barrier can 
strongly influence the transport of discrete solitons. It 
is therefore of interest to investigate the conditions un- 
der which the discrete soliton ratchets can exist. The 
present paper is just devoted to this investigation. More 
precisely, we study the ratchet dynamics induced by tem- 
porarily asymmetric forces of zero mean on topological 
solitons (kinks and antikinks) of the discrete sine-Gordon 
(DSG) system. This equation models a number of physi- 
cal systems such as arrays of Josephson junctions, crystal 
dislocations or charge-density waves (see [23,24]). In par- 
ticular, we investigate the conditions for the occurrence of 
soliton ratchets and study the dependence of the average 
soliton velocity on the system parameters. A comparison 
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with the results derived for continuous soliton ratchets in 
Ref. [10] is also provided. 

From our study it emerges that discrete soliton ratch- 
ets are much more complicated than their continuum 
counterpart. In particular, the mean velocity of the kink 
in most cases appears to be a piece- wise function of the 
parameters which resembles a devil's staircase. We find 
that kink transport becomes very effective on the cor- 
responding orbits (limit cycles) which are phase locked 
to the external driver. Transport is possible also in the 
presence of more complicated dynamics such as chaotic 
and intermittency orbits, especially when the system be- 
comes very discrete (this is achieved when the coupling 
constant is very small). In these cases, however, the kink 
transport is not very efficient since the drift velocity is 
rather small. Except for the very discrete case, domi- 
nated by the pinning of the kink to the lattice, we find 
that the internal mode mechanism remains valid also in 
the DSG ratchet in all cases in which the transport is 
observed. 

We remark that the soliton ratchets induced by tem- 
porarily asymmetric fields may be an effective way to con- 
trol the transport properties of a large variety of contin- 
uous and discrete systems. From the experimental point 
of view, indeed, it is much easier to produce ratchets by 
means of asymmetric fields than by breaking the internal 
spatial symmetry of the system [25] , since in the former 
case no structural changes of the system are required. 

The paper is organized as follows. In Section II, we 
present the model and derive the necessary condition for 
the directed kink motion in terms of a simple symmetry 
analysis, which is based on a point particle description of 
the kink dynamics. In this section we also describe the 
desymmetrization mechanism and confirm the results by 
means of numerical simulations. In the next section, we 
study the dependence of discrete kink ratchets on the sys- 
tem parameters. In particular, we investigate the depen- 
dence of the mean velocity of the kink on the amplitude, 
phase, and frequency of the AC driver as well as on the 
damping and the coupling constant. Moreover, the valid- 
ity of the internal mode mechanism in the discrete case is 
discussed. In Section IV, we consider the soliton ratchet 
in a finite lattice and discuss possible applications of the 
phenomenon to arrays of small Josephson junctions. Fi- 
nally, in Section V, we summarize the main conclusions 
of the paper. 

II. THE MODEL 

The AC-driven and damped discrete sine-Gordon 
(DSG) equation is introduced in a dimcnsionlcss form 
as follows 

u n — k Au n + smu n + aii n + E{t) =0, n = 1, 2, . . . N. 

(1) 



Here u n is the displacement of the nth particle from its 
equilibrium position, Au n = u n+ \ — 2u n + is the 

discrete Laplacian, k is the coupling constant measuring 
the discreteness of the lattice, a is the damping coefficient 
and E(t) is an external driving field. In the following we 
assume E(t) to be of the form 

E(t) = Ex cos(ui) + E 2 cos(mujt + 6) , (2) 

where m is an integer even number. Notice that the su- 
perposition of two harmonics makes the periodic force to 
be asymmetric in time for almost all values of 6, a fea- 
ture which can be used to break the temporal symmetry 
of the system (see below). In this context, it is of in- 
terest to investigate the condition under which a driving 
force of zero mean of type (2) can induce kink's unidi- 
rectional motion similar to the one observed in the con- 
tinuum SG case [10]. In this regard, we remark that 
in the lower approximation, a discrete kink of the form 
u n (t) = 4arctan{cxp[n — X(t)]}, can be viewed as a sin- 
gle particle [24] and its dynamics is described in terms 
of collective coordinates: the center of mass X(t) and 
the kink velocity X{t). In this approach, the effective 
point-particle equation of the motion becomes 

X + aX + V^ N (X)+E(t) = , (3) 

where V PN {X) = V PN (X + l), V PN (X) ~ sin27rX is the 
Peierls-Nabarro (PN) potential accounting for the dis- 
creteness of the lattice and E(t) <~ E{t) is the effective 
driving field of the kink [we assume E(t) to be propor- 
tional to E(t)]. An important parameter of the problem 
is the frequency of kink oscillations in the bottom of the 
PN potential (the PN frequency), which can be written 
as [30] 

ujpn = \J 27ra K 3 / 2 exp (-7T 2 \/k/2) , a ~ 307r . (4) 

Within this approximation, the unidirectional motion of 
the kink corresponds to a limit cycle of Eq. (3) , which is 
phase locked to the frequency of the external driver. On 
this orbit, the average kink velocity is expressed as 

<«> = <*(*)> = y ' £ ■ (5) 

with k and I being integer numbers. Notice that in this 
resonant regime, the kink travels k sites during the time 
IT = 2irl/cu, so that, except for a shift in space, its profile 
is completely reproduced after this time interval (in the 
pendulum analogy, this orbit corresponds to k full rota- 
tions of the pendulum during I periods of the external 
drive). In the following, we will refer to the phase locked 
dynamics also as to resonances. 

A. Symmetry properties and conditions for 
transport 

In analogy with the continuous SG case [10,11], one can 
expect that the directed kink motion arises when all the 
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symmetries of Eq. (3), which relate kink solutions with 
opposite velocities, are broken. Qualitative conditions for 
the occurrence of this directed motion can be obtained 
from the analysis of the symmetry properties of Eq. (3). 
In this approach, the many-particle problem is reduced to 
the one-particle ratchet studied before [3,4] [we neglect 
oscillations of the discrete kink profile (to be discussed 
later) which also contribute to the phenomenon]. Our 
analysis is based on the simple observation that the sign 
of the soliton velocity X(t) can be changed by means of 
the following symmetry operations: 



T> x :t^t + T/2, X -> -X, 

D r : X -> X + X , t -> -t + 2t , 



(6) 
(7) 



where t)x denotes a shift in time followed by a reflection 
is space and E>r is a shift in space followed by a reflection 
in time (here to is a constant and Xq is either an integer 
or a half integer) . Notice that Eq. (3) is always invariant 
under the symmetry T>x provided the external driver 
satisfies 



E{t) = -E(t + T/2) 



(8) 



Notice that Vpn{— X) = Vpn(X) is always satisfied be- 
cause the sine-Gordon potential is symmetric. In the 
zero damping limit (a — > 0), Eq. (3) becomes invariant 
also under the symmetry T)q- with Xq — 1, provided the 
external driver satisfies the condition 



E(t + t )=E(-t + t Q ) 



(9) 



with to being a constant, which depends on the shape of 
E(t). In the overdamped limit (a — * oo), Eq. (3) can be 
rewritten as X + V PN (X) +E(t) = from which one can 
see that it becomes invariant under the symmetry D7- 
with Xq = 1/2, provided the external driver satisfies the 
condition 



E(t + t ) = -E(-t + t ) 



(10) 



It should be remarked here that V PN {X) 



-v PN (x 



1/2) is always satisfied since for the DSG equation we 
have that V PN (X) ~ sin27rJT). From the above symme- 
try properties it follows that one can break all symme- 
tries relating orbits with opposite velocities by properly 
choosing the driving force E(t). Thus, for the general 
case a ^ 0, we have only the symmetry T>x for Eq. (3). 
This symmetry can be broken by choosing any function 
E(t) which violates Eq. (8). In the zero damping limit 
(q — > 0), we have, besides D^, also the symmetry D7- 
with Xq = 1. In this case, one must choose a driving field 
E(t) which violates both Eq. (8) and Eq. (9). Similarly, 
in the overdamped limit (a — » 00), a function E(t) which 
violates both Eq. (8) and Eq. (10) should be chosen. 

From these considerations it follows that a simple si- 
nusoidal driver cannot support the kink transport in the 
lattice, since Eqs. (6)-(10) in this case are always satis- 
fied. For a biharmonic driver of the type (2), however, 



we have that Eq. (8) is violated for any 9 (if m is even), 
so that the kink transport should become possible. 
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FIG. 1. Contour plot of the temporal evolution of the par- 
ticle velocities u n (t) for (a) E\ = 0.3, E2 = 0.15 and (b) 
Ei — 0.45, E2 — 0. Other parameters are: k = 1, a = 0.05, 
uj = 0.1, 9 = 2, D = 0.002, N = 500. Periodic boundary 
conditions have been applied. 

From these arguments it is also clear that an external 
periodic driver of zero mean, which consists of the super- 
position of only the first two harmonics, is the simplest 
driver that can be used to induce the soliton ratchets in 
the DSG system. 



B. Numerical study of transport vs symmetries 

To verify the validity of the previous analysis and to 
check the desymmetrization of the orbits as a function of 
the driver parameters Ei,E%, and 8, we recourse to direct 
numerical integration of Eqs. (1). In order to be sure 
that the system explores the whole phase space and that 
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the phenomenon does not depend on initial conditions, 
we perform the first step simulations in the presence of 
white noise. 



investigations (see below and Sec. IV) show that this is 
true in almost all cases, except for very narrow intervals 
in the parameter space where two attractors can coexist. 
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FIG. 2. Dependence of the average kink velocity on the 
phase difference 8 for E\ — E2 — 0.2, a — 0.05 (a), 
a = 0.2 (b) and a — 0.5 (c). Other parameters are: K = 1, 
uj — 0.35. Corresponding pairs of rotation numbers (k, I) are 
given nearby the most pronounced resonances. Dashed lines 
mark the coordinate axes. 

White noise has been included by adding in the r.h.s. 
of Eqs. (1) a stochastic term £ n (t) of zero mean, (£ n (i)) — 
0, and with the autocorrelation function (£, n (t)£, m (t')) = 
2aD5 mn S(t — t'). The resulting Langevin equations have 
been integrated numerically by using the fourth-order 
Runge-Kutta method, adopting either free ends or peri- 
odic boundary conditions: u„ + jv(rj) = u n (t)±2ir (positive 
and negative signs refer to kinks and antikinks, respec- 
tively). 

Figure 1 illustrates the dynamics of a kink of the 
damped DSG equation under the influence of noise, 
driven by a biharmonic driver (2) with m = 2 [panel 
(a)], and by a single harmonic driver [panel (b)]. We see 
that while the single harmonic driver is unable to produce 
directed motion, the biharmonic driver is quite effective 
to produce the kink transport. From this figure, one can 
also observe that there is only one attractor correspond- 
ing to the unidirectional motion of the kink. Further 
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FIG. 3. Velocities of two kinks, running into the opposite 
directions at a = 0.05, u — 0.1, E\ = 0.21, 6 — as a function 
of E2/E1. Dashed line denotes that the given solution is no 
longer stable. 

This fact makes possible to investigate Eqs. (1) in the 
absence of noise without averaging over the initial con- 
ditions. In Fig. 2 we plot the time average kink velocity 
(v) as a function of the phase difference 9 in the case 
of a biharmonic driver with m = 2. In contrast to the 
continuous case for which (v) was shown to have a sinu- 
soidal dependence on 9 [10,12], in the discrete case, we 
find a complicated piece- wise dependence (v){9), which 
resembles a sinusoidal function only slightly. This is due 
to the fact that in the most of cases the dynamics is 
phase locked to the external driver and the kink velocity 
is given by Eq. (5). 

Notice that for a weak damping [see panel (a) of Fig. 2] 
several resonances (k = I = ±1, k = ±2,1 = 1) are 
clearly visible, while for a strong damping [panel (c)], 
the kink becomes pinned to the lattice for large intervals 
of 9 (depinning of the kink would require stronger fields). 
From Figs. 2(a-c) a qualitative understanding of the sym- 
metries to break, in order to achieve unidirectional mo- 
tion, can be obtained. In particular, one can see that 
close to the underdamped limit (a = 0.05), the mean ve- 
locity becomes zero in the intervals 9 6 [—2.26,-1.81], 
9 e [0.80,1.06], and 9 e [1.11,1.33] [see Fig. 2(a)]. For 
larger values of the damping constant [see Fig. 2(b)], we 
have (v) = in the intervals 9 G [—2.66,-1.88] and 
9 £ [0.48, 1.26], while for a = 0.5, the kink stays pinned 
for 9 e [-2.5,-0.95] and 9 S [0.64,2.19]. By increasing 
a, the extremal values of (v) shift to or ±7r with the 
growth of a. Moreover, the length of the intervals in 9 
where (v) — increases, while the average kink velocity 
decreases as a increases. This behavior is an obvious con- 
sequence of the slowing down effect of the damping on the 
kink motion. In the overdamped limit a — ► 00, the in- 
tervals become centered around the points 9 = ±7r/2 for 
which E{t) = E x cosujt + E 2 cos(2ujt + ir/2) = -E(w/ui- 
t) and thus the symmetry (6) is satisfied. Instead, in the 
underdamped case (a — > 0), the centers of these inter- 
vals gradually shift to the positions 9 — 0, ±ir for which 
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E(t) = E(—t) and the symmetry (7) is satisfied. 

We recall that in the continuum case [10], the func- 
tional dependencies obtained in the underdamped and 
overdamped limits were found to be (v) ~ sin# and 
(v) <~ cos 6*, respectively. Thus, the transition from one 
limit to another in the discrete case is similar to the con- 
tinuous case in the sense that the directed soliton motion 
disappears in correspondence with values of 9 for which 
the respective symmetry is restored. 
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FIG. 4. Dependence of the average kink velocity on the 
amplitude of the driver Ei — E2 for 9 = —0.5 (decreasing 
dependence) and 9 — 2 (increasing dependence). Other pa- 
rameters are: a = 0.1, k = 1, ui = 0.35. Numbers in brackets 
show the respective pair of rotation numbers (k,l). Insets 
show more details for case 9 = —0.5. 

In the discrete case, however, the soliton velocity is 
zero not only for those values of 9 which restore the 
symmetries (6) and (7), but also for some finite inter- 
val around these values. 

In the case of odd values of to, the condition (8) is sat- 
isfied for any value of 9 if E 2 ^ 0. We have indeed that 
E(t + T/2) = E x cos(wi + n) + E 2 cos(mujt + tott + 9) = 
—E(t) and therefore the symmetry T)x is always valid. 
This observation implies that there should be no kink 
transport in the system, a result which is indeed con- 
firmed by direct numerical simulations of the full system 
with to = 3, both in the presence and in the absence 
of noise. Also, in analogy with the continuous case, we 
find that for the same parameter values the discrete an- 
tikink ratchets always move in the direction opposite to 



the kink motion. 

The emergence of the directed kink motion can also 
be seen as a consequence of the desymmetrization of the 
basins of attraction of the two limit cycles corresponding 
to kinks moving with opposite velocities (for single parti- 
cle ratchets, the desymmetrization of the orbits has been 
shown in Rcfs. [3,31]). 
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FIG. 5. Coordinate of the kink center X c as a func- 
tion of time for different values of the driving amplitude 
E 1 = E 2 = 0.129 (a); 0.12981 (b,c) and 0.1309(d,e). In panels 
(c) and (c), X c (t) has been plotted after each oscillation pe- 
riod T (see text for details). Other parameters are: a — 0.1, 
w = 0.35, k = 1, 9 = -0.5. 

As soon as the symmetry is broken (by switching on 
the field E 2 ), one of the basins of attraction begins to 
shrink and eventually disappears as E 2 increases. As a 
result, only the limit cycles which correspond to motion 
in one direction survive, as one can see from Fig. 3. In 
particular, for the parameters given in this figure, we find 
that the attractor corresponding to the limit cycle with 
rotation numbers k = 1,1 = 1 disappears at E 2 /Ei — 
0.0118, i.e., already for a rather weak asymmetry of the 
field (E 2 < E^. 



III. DEPENDENCE OF THE PHENOMENON ON 
SYSTEM PARAMETERS 

In this section, we investigate the dependence of kink 
transport on the system parameters and the internal 
mode mechanism. 



A. Dependence on the driving amplitude 

For continuous ratchets it was shown previously [10] 
that the average kink velocity is proportional to E\E 2 , 
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so that, provided the respective symmetries are broken, 
the directed motion can occur for arbitrary small values 
of the driving amplitudes. For the DSG equation, the 
dependence of the average kink velocity on the driver 
amplitudes is depicted in Fig. 4 from which one can see 
that the kink velocity is not a smooth monotonic function 
of £a, 2, but a piecewise function with plateaux of differ- 
ent lengths, resembling a "devil's staircase" [the same as 
in Fig. 2(a) but monotonic]. The plateau values of the 
kink velocities are given by Eq. (5) and correspond to 
dynamical regimes, which are limit cycles with rotation 
numbers (k, I) phase-locked to the driver. Notice that the 
largest resonant step ((v) = 0.0557) is achieved at the 
rotation numbers corresponding to the main resonance 
(k = ±1,1 = 1). One can easily observe the smaller 
resonant steps with k — ±2,1 = 1 and k = ±1,1 = 2, 
for which the kink velocities equal (v) — 0.1114 and 
(v) = 0.02785, respectively. Higher order resonances can 
also be identified, but they are not well visible in the 
figure. 

The main feature of the discrete case is the existence 
of a threshold in the driver amplitude (depinning thresh- 
old), above which the transport can occur and below 
which the kink is pinned to a lattice site, oscillating 
around its center of mass: 



In Fig. 5, the kink position X c is depicted as a function 
of time for different values of the driving amplitude. We 
see that when the driving amplitude is below threshold 
[panel (a)] , the kink performs periodic oscillations around 
its center of mass, while above the threshold [panels (b)- 
(e)], the standing kink becomes unstable and the directed 
motion starts. By further increasing the driving ampli- 
tude above the threshold, the kink motion becomes either 
phase locked to the external driver or chaotic with an in- 
termittent behavior [see panels (b) and (c)]. Panels (d) 
and (e) of Fig. 5 show the phase-locked dynamics with 
rotation numbers k = 1,1 = 5. Notice that in this case, 
the kink travels over a fixed number of sites during each 
period T in agreement with Eq. 5. By further increasing 
E\, the the kink dynamics starts to switch between pe- 
riodic (or quasiperiodic) regime and intermittency. Sim- 
ilar transitions in the case of a single harmonic driver 
nave i ilso been reported in Refs. [23,28]. In the intermit- 
tency regime, the kink dynamics switches in an unpre- 
dictable manner between two attracting limit cycles [see 
panels (b) and (c) of Fig. 5]. In this case, the kink dwells 
some time interval on a certain site, pinned by the corre- 
sponding minimum of the PN barrier, before jumping to 
the next site. From panel (b) one can also see that the 
dwelling time is not rationally related to the driving pe- 
riod T and changes randomly from jump to jump. Thus 
in the intermittent regime the kink does not travel the 
same number of sites during an integer number of periods 
of E(t) and the global dynamics is chaotic. 



The same behavior is seen in panel (c) of Fig. 5 for a 
larger time scale. Note that the intermittency behavior 
occurs not only around the depinning threshold, but also 
for larger values of the driving amplitude. Also notice, 
from the inset (b) of Fig. 4, that the depinning threshold 
has an hysteretic behavior. In this regard, we remark 
that the numerical investigation has been performed by 
increasing the driving amplitude in small steps, taking 
the final state for a given step as an initial condition 
of the next one. When the pinned state loses stability, 
the system finds itself on a chaotic attractor, which cor- 
responds to the directed motion. By increasing further 
Ei we obtain larger values of the kink velocity, while if 
we move backwards, we observe that two attractors can 
coexist - one of them corresponding to a pinned kink 
oscillating around its center of mass, and the other to 
a kink performing ratchet dynamics. If Ei is further 
decreasing, the moving state loses stability and the sys- 
tem jumps back to the pinned state. The width of the 
hysteresis (interval between two bifurcations) appears to 
be rather small, i.e., 0.129065 < E 1 < 0.12979. Simi- 
lar hysteretic phenomena also appear at larger values of 
the driver amplitudes. For example, one can find that 
in the interval 0.1660 < E\ < 0.1678, two limit cycles 
with k = 1, 1 = 2 and k = 1, 1 = 1 coexist. A more de- 
tailed numerical investigation of the (v)(Ei) dependence 
shows that the observed devil's staircase is incomplete 
and there are gaps inside it, which depend on the value 
of 9. As the (v)(Ei) dependence becomes steeper, less 
phase locked states are found. In particular for the case 
9 = 2, we see that in the transition from the k = 1, 1 = 2 
to the k — 1,1 — 1 state, all intermediate rational values 
of k and I are missing. 

B. Dependence on the driving frequency 

The dependence of the mean kink velocity on the driver 
frequency was investigated in Ref. [10] for the continu- 
ous SG model. In this case, it was shown that the aver- 
age velocity depends on the driving frequency as (v) ~ 
sin [9 — 9 (uj; a)] /uj 3 , where tan6>o(w; a) = [a/(2ui)] ■ [3 + 
(a/uj) 2 ]. This fact implies that max e / e [ .2 7r [( w ) decays 
with the growth of lo. For a fixed value of 9, however, this 
dependence is defined by the mutual relations between ui 
and a, so that (v) can experience oscillations and sign 
reversals before tending to zero either in a decreasing or 
in increasing way. Similar behavior is expected also for 
the DSG equation, although in this case, the problem is 
complicated by the presence of the Peierls-Nabarro (PN) 
frequency lupn in Eq. (3). In analogy with Fig. 2, one 
can expect the dependence (v)(lo) to be also a piecewise 
function with a "devil's staircase" character. 

In Fig. 6, we depict the kink velocity, normalized to 
the driving frequency ((v)T — 2-itv/uj) as a function of oj 
for two values of the damping constant (notice that for 
phase locked dynamics, (v)T coincides with the ratio of 
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rotation numbers k/l). One can see that the dynamics 
is characterized by a series of resonances, the most pro- 
nounced one being at the main frequency (k = 1,1 = 1). 
We also observe that the resonances with k > I are less 
pronounced and the subharmonic resonances with k < I 
are practically not visible in the figure. The case of small 
damping and larger driver amplitude (depicted in the 
figure by circles) is characterized by an almost mono- 
tonic decay, while in the opposite case (larger damping 
and smaller driving amplitude), there are wide pinned 
regions and some peak around the value of lopn [notice 
from Eq. (4) that lj pn ~ 0.17 for n = 1]. 



while inset (c) shows details of the (v) dependence in the 
neighborhood of the second resonance. The coupling of 
the kink motion with linear waves occurs in the interval 
0.447 < lo < 0.517. At the beginning of this interval, 
the kink displays a chaotic tail as shown in Fig. 6(a). As 
frequency increases, the chaotic tail becomes more and 
more regular. This is shown in inset (b) of the figure, 
from which one can see that at lo — 0.507 the oscillating 
tail is almost monochromatic. At lj = 0.517 the oscil- 
lating tail turns into a localized oscillating mode which 
decays at infinity. A further increase of the driving fre- 
quency causes the decay of the width of the mode. 
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FIG. 6. Dependence of the normalized kink velocity (v) T 
on the driving frequency for Ex = E 2 = 0.1, = 2, a = 0.05 
(o) and Ex = E 2 = 0.02, 9 = 1.5, a = 0.15 (x). In both 
cases k = 1. The insets (a) and (b) show the kink profile 
for the first set of parameters when lj = 0.45 and lj = 0.507, 
respectively. Inset (c) shows more detailed behavior around 
lj = 0.5. 

In both the cases we find that, except for the reso- 
nances lj i=a lj l and 2lj w lol, where is the frequency 
of linear waves 



;£(«) 



2k(1 — cosq) + 1, 



(12) 



the kink transport becomes effective mainly at low fre- 
quencies lo -C ljl and disappears for driving frequencies 
lo > lol- Also, from Fig. 6 one can see that the kinks 
become pinned to the lattice for values of lo significantly 
smaller than ljl- At the resonances and 2lo ~ lol, 

we observe that the kink dynamics become coupled to lin- 
ear waves in the system (plasmons in the case of array of 
Josephson junctions). Insets (a) and (b) of Fig. 6 show 
the profiles of the kink solution in these resonant cases, 
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FIG. 7. Time evolution of the position of the kink cen- 
ter for: (a) k = 0.5; Ex = E 2 = 0.21 (curve 1), 0.19 
(curve 2), 0.17 (curve 3), 0.11 (curve 4) and (b) k = 0.25; 
Ex = E 2 = 0.14 (curve 1), 0.24 (curve 2), 0.22 (curve 3). 
Other parameters are: 6 = 1.5, a = 0.1, lj = 0.35. The 
upper inset of panel (a) corresponds to details of the case 
Ex = E 2 — 0.19 and the lower inset corresponds to the case 
Ex = E 2 = 0.11. The inset of panel (b) corresponds to the 
case Ex = E 2 = 0.22. In all figures, the data have been 
plotted with the interval T = 2-k/uj. 

A similar scenario occurs also at lo rs lol, but in this 
case the kink velocity is much smaller and practically not 
visible in Fig. 6. Notice that in the frequency interval 
0.895 < lo < 1.05, the kink becomes again coupled to 
linear modes and displays a chaotic tail, which becomes 



7 



more regular as u) further increases. At the beginning 
of this window, the dynamics in the tails are strongly 
chaotic and they are accompanied by the formation of 
large- amplitude localized excitations (breathers). 

By reducing the driving amplitude, the frequency win- 
dows, in which the coupling with linear waves occurs, 
decrease. We find that for E\ = £2 = 0.04, the first 
coupling window occurs at 0.480 < u < 0.495 and the 
second one at 0.955 < u> < 0.994. In this case, no chaotic 
dressing of the kink is observed [the coupling occurs with 
very few (or single) linear modes and the kink looks very 
similar to that in Fig. 6(b)]. For the driving amplitude 
Ei = £2 = 0.01, the coupling with linear waves does not 
take place at all, and around the resonant frequencies 
u>l and ujl/2, the kink remains pinned to the lattice. As 
the driving amplitude increases, beyond a certain thresh- 
old chaotic oscillations can completely destroy the kink. 
From these results, we conclude that the coupling of the 
kink dynamics with linear waves depends very much on 
the amplitude of the external field and this effect is non- 
linear in E\. Near the resonances with linear modes, 
the kink dynamics become more chaotic (diffusive), re- 
sembling the one of a Brownian particle (the kink makes 
many random jumps backward and forward and the uni- 
directional motion can be seen only at a large time scale) . 

C. Dependence on the coupling constant 

In this section, we investigate the dependence of the 
ratchet phenomenon on the interaction constant k and 
the conditions to sustain mobile kinks in the system. As 
well known [32], the discreteness usually prevents free 
topological soliton propagation. In our model, the dis- 
creteness of the system is characterized by the interaction 
constant k, with k — > 00 corresponding to the continuum 
limit. In Fig. 7, we depict the temporal evolution of the 
kink center of mass X c for two values of the coupling con- 
stant k = 0.5 [panel (a)] and k = 0.25 [panel (b)]. From 
this figure, we conclude that the unidirectional motion 
exists also for small values of k with dynamics which is 
chaotic rather than phase locked (see insets of the figure) . 

Chaotic motion, leading to transport is found to be 
either of the intermittent type [curve 4 of panel (a)] or 
diffusive type [curves 1-3 of panel (a)]. The intermittent 
regime is more often observed for k > 1 and the behav- 
ior is the same as described in the previous subsection, 
i.e., the kink oscillates in the minima of the PN poten- 
tial until it jumps to the adjacent site [see lower inset 
of panel (a)]. The motion is characterized by the fact 
that the jumps of the kink occur randomly in time, but 
always in the same direction. For k < 1 the diffusive mo- 
tion is the most typical scenario. In this case, the kink 
jumps randomly forward and backward [see upper inset 
of panel (a) and inset of panel (b)], but on average the 
motion remains unidirectional. A decrease of the cou- 
pling constant makes the dynamics even more diffusive 



in the sense that the forward and backward jumps be- 
come of larger amplitude with a consequent decrease of 
the average velocity. 

In the left panel of Fig. 8 the diagram of possible dy- 
namical regimes in the plane (k,E\) is shown. We ob- 
serve that by changing the coupling constant k, one can 
pass from regular (phase locked) dynamics to chaos. The 
regular dynamics are dominant for n > 1 and when k 
is decreases, the windows of chaotic motion appear. By 
further decreasing k, the number of chaotic windows in- 
creases and the chaotic motion becomes dominant. At 
k = 0.25, almost no regular dynamics exist. Notice that 
the diagram in Fig. 8 does not show the complete pic- 
ture and the details of transitions to different dynamical 
regimes because in order to limit the computation time, 
we have used in the numerical calculations the step of 
0.01 in the amplitude of the driving field E\ — E2 and a 
coarse step of 0.25 in k. We find that for driving ampli- 
tudes Ei = E 2 > 0.3, the chaotic dynamics of the whole 
lattice destroy kink solutions. 
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FIG. 8. Left panel. Diagram that demonstrates main 
regimes of the kink ratchet motion for different values of the 
coupling constant k and the driving amplitude E\ — E 2 (other 
parameters are: 9 — 1.5, a — 0.1, and u) = 0.35). In the dia- 
gram o stands for the regime when kink is pinned ((v) = 0), 
© stands for the periodic or quasiperiodic ratchet motion, x 
stands for the chaotic regime, and ♦ corresponds to the case 
when kink solutions do not exist. 

Right panel. Dependence of the average kink velocity on the 
driving amplitude E\ = E2 for different values of the cou- 
pling constant k = 0.5(A), k = 0.75 (x), k = 1 (o), k = 1.5 
(♦), and k = 2 (o). Other parameters are as in the left 
panel. The insets show dependence of the average velocity 
on the phase difference 0. Inset (a) corresponds to k = 0.75, 
Ei = E 2 = 0.19 (o) and Ei = E 2 = 0.25 (x), inset (b) corre- 
sponds to k — 0.5, Ei — E 2 — 0.22 (©). Solid lines are drawn 
as a guide to an eye. 

In the right panel of Fig. 8, we show the dependence of 
the average kink velocity on the driving amplitude E\ = 
E2 for different values of the coupling constant k. We 
observe that already for n — 2, the curve is very smooth 
(except for the small resonant steps k = 3, 1 = 1 and k = 
3,1 = 2) and the behavior becomes very similar to that 
reported for the continuous limit ((v) ~ EiE 2 ) [10]. By 
increasing the driving amplitudes £1,2, the dependence 
becomes non-monotonic, a fact which was also observed 
in the continuous SG case and ascribed to the interaction 
of the kink with internal oscillation modes [10]. On the 
other hand, decrease of the coupling constant makes the 
phase locking steps more pronounced. For n = 1.5, the 
dependence (v)(Ei) almost coincides with the one for 
k = 2 (for larger amplitudes, however, the three resonant 
steps k = I = 1, k = 3, 1 = 2, and k = 2, 1 = 1 become 
much more visible). 

Further decrease of the coupling constant makes the 
dependence even less monotonic: for n = 1 (shown by o 
in the figure), the non-monotonicity is quite weak, with 
a small interval 0.145 < E\ < 0.154 in which the kink ve- 
locity drops back to zero. For E\ > 0.2, the dependence 
is also non-monotonic because of the interaction of the 
kink with vibrational modes localized on it. For k = 0.5 
and k = 0.75, the average kink velocity becomes strongly 
non-monotonic. 

Notice from the right panel of Fig. 8 that several re- 
versals of the kink motion along the E\ axis occur. In 
this regard, we remark that the sign of (v) depends on 
the relations between u), a, and uopn- Decrease of n can 
change lopn significantly [see Eq. (4)], thus effecting the 
sign of the velocity. The reversals of (v) as a function of 
Ei, however, are not fully explained by these arguments 
and we believe that, in analogy with the continuous case, 
the coupling of the kink with small amplitude waves plays 
also an important role. 

By comparing the left and the right panels of Fig. 8, 
one can also see that the non-monotonic jumps occur 
in correspondence with the transitions from chaotic to 
regular (phase locked) regimes. In order to get a better 



understanding of this behavior, we have plotted in insets 
(a) and (b) the mean velocity as a function of the phase 
difference 9. The inset (a) refers to the case k = 0.75 for 
two different values of the driver amplitude. We see that 
while for E\ = 0.19 the dependence (v)(6) is piecewise, 
but somewhat similar to Fig. 2 (i.e., it has one minimum 
and one maximum), the dependence becomes more irreg- 
ular as Ei increases, with small islands of phase locked 
regimes with opposite velocities. A similar phenomenon 
is observed, when n decreases to k = 0.5 [inset (b)]. 
This behavior is linked to the well known phenomenon of 
crises of attractors, i.e., to sudden appearance (or disap- 
pearance) of an attractor as a system parameter varies 
(see [33]). By further decreasing the coupling constant, 
more and more attractor crises are found. 

From the insets of Fig. 8, we also see that, in contrast 
with the case n = 1 in Fig. 2, the dependence (v) — (v)(9) 
loses any resemblance with a sinusoidal function for small 
values of n. From Fig. 8, it is also evident that the de- 
pinning amplitude of the kink motion increases when k 
decreases. For very small k, e.g., n = 0.0625 (not shown 
in Fig. 8), no directed motion is found - the kink is ei- 
ther pinned by the lattice or destroyed by chaos when the 
driving amplitude becomes large enough. This behavior 
is not surprising, since we know that high discreteness 
normally prevents solitons from free propagation. 




FIG. 9. Evolution of the kink profile u n (t). Upper panel 
shows initial kink profile (+) and kink profile after integration 
time 7T (o), shifted two sites backwards. Parameters of the 
model are: Ei = E 2 = 0.135, 6 = 1.5, a = 0.1, to = 0.35, and 
k = 1. 

We remark that the attractor crises leading to reversal 
phenomena are out of the range of validity of the single 
particle approximation since they involve many degrees 
of freedom. In particular, in the limit of weak coupling, 
the point particle approximation (3) as well as the sym- 
metry analysis discussed in Sec. II become not valid. 
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D. Internal mode mechanism and dependence on the 
damping parameter 

In the case of continuous soliton ratchets, it was shown 
(see Refs. [9,10,19,20,18]) that a contribution to the 
ratchet phenomenon comes also from the internal oscil- 
lation of the kink via the internal mode mechanism. We 
expect this effect to be true also in the present case, both 
for large values (i.e., close to the continuum limit) and 
intermediate and small values of the coupling constant k 
(for very small values of k, however, the kink becomes 
pinned to the lattice and the ratchet phenomenon disap- 
pears as discussed above). 

In the following we investigate the internal mode mech- 
anism by fixing k = 1 and performing direct numerical 
simulations of Eqs. (1). In particular, we show the exis- 
tence of a local oscillation on the kink profile which is, 
perfectly synchronized (phase locked) with the kink mo- 
tion. We have found that when the dynamics of the kink 
center of mass is phase locked to the resonance (k,l), 
the internal mode oscillation is also locked to the same 
resonance. This is shown in Fig. 9 for the discrete kink 
ratchet, which is phase locked to the external driver with 
rotation numbers k = 2, I = 7. In the panel (a) of 
this figure, the initial and final kink profiles are depicted, 
from which one can conclude that these profiles perfectly 
coincide after I = 7 periods of the external driver (notice 
that the final configuration has been shifted by two sites 
backwards in order to demonstrate full coincidence). 
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FIG. 10. Same as in Fig. 9. All parameters are the same 
except Ei = E2 = 0.17. In panel (a), the initial kink config- 
uration is shown by (+), the configuration after time T (o) 
and after time 2T (o), shifted by one site backwards. Panel 
(b) shows three-dimensional dynamical picture for the kink 
evolution in the time interval [0, 4T]. 

The phase locking of the motion to the external driver 



is also observed from the three-dimensional plot of u n (t) 
in panel (b). Similar results arc obtained for different 
values of the driving amplitudes, as one can see from 
Fig. 10. In this case, the kink is locked on the resonance 
with rotation numbers (1,2) (notice that the kink repro- 
duces itself completely after the time 2T). In panel (b), 
the complete dynamics in the time interval [0, 4T] is also 
shown. The fact that the center of mass motion and 
the oscillations on the kink profile (internal mode) are 
perfectly synchronized suggests the existence of coupling 
between the internal and translational modes similarly 
as for the continuous case. 

The existence of the internal mode mechanism in the 
discrete case is also supported by the influence of the 
damping constant on the phenomenon. In the continuous 
case, it was shown that the coupling between the trans- 
lational and internal modes decreases with the damping 
and for the case of asymmetric potentials with symmet- 
ric drivers, it was proved that in the limit a — > 0, the 
coupling completely disappears [18]. In the case of bi- 
harmonic asymmetric forces, the nonlinearity induces an 
effective bias component, which gives rise to a point par- 
ticle contribution also in the absence of damping. 
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FIG. 11. Existence diagram of the main resonance 
(k = I = 1) on the plane (a, Ei = E 2 ) for ui = 0.35, 9 = 2, 
and k = 1. The inset shows the dependence of the width of 
resonance AE on damping coefficient a. When the resonance 
consists of several islands, AE is computed as a sum of the 
widths of individual islands. 

The contribution of the internal mode to the ratchet 
dynamics, however, should be sensitive to the damping, 
if the coupling is controlled by the damping and the uni- 
directional motion should become less effective for small 
values of the damping constant. On the other hand, when 
the damping in the system becomes too large, the dy- 
namics are strongly reduced (or stopped) , so that a non- 
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monotonous behavior of transport with a is expected and 
an optimal value of the damping which maximizes the 
transport should exist. 

In the following numerical study we have taken the 
width of resonance steps produced by the ratchet dy- 
namics (see Fig. 4) as a measure of the efficiency of the 
transport in the system (in the experimental setting con- 
sidered in Rcf. [14], these resonance widths correspond to 
the voltage steps in the IV characteristic of a Josephson 
junction). 

Figure 11 illustrates the width of the main resonance 
(1,1) as a function of a (similar behavior is observed 
also for other resonances). It follows from this figure 
that the resonant steps are reduced when a is reduced 
and tend to disappear as a — > (for small a numerical 
calculations become very difficult due to longer transient 
times to settle on the attractor). The same behavior is 
observed as a increases beyond the value a w 0.35. From 
this figure we conclude that the effect becomes maximum 
around an intermediate value of damping with a m 0.2. 



the nth junction. The discreteness parameter equals 
k = y $o/(27r/ c L), where $0 is the magnetic flux quan- 
tum, L is inductance of an elementary cell, and I c is the 
critical current of an individual junction. The dimcn- 
sionless dissipation parameter is then a = $ /(2irI c R), 
where R is the resistance of an individual junction, and 
the time is normaliz ed to the inv erse Josephson plasma 
frequency 1/ljq = \JC<&q/{2ttI c ). In these systems, the 
topological solitary waves have the physical meaning of 
trapped magnetic flux quanta (fluxons) and the voltage 
drop V in the array is defined as 

V = ^J2 , lim 7 (13) 

The experiments with annular JJAs have been per- 
formed for typical lengths TV ~ 8 — 30 (see Refs. [34-36]). 
In the following, we consider the case of an array with 
N = 10 junctions subjected to periodic boundary con- 
ditions (annular array). From the previous analysis we 
expect that the ratchet dynamics of the fluxon give rise 
to nonzero voltage drops in the IV characteristic of the 
array. 
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FIG. 12. Dependence £1 = E^V) for n = 0.25 (o) and 
k = 1 (x). Other parameters are: a — 0.1, 9 — 1.5, and 
u) = 0.35. The array consists of N = 10 junctions. Dashed 
line for the case k = 0.25 has been used as a guide for an eye. 

The existence of an optimal coupling between the kink 
phase looked dynamics and its internal oscillation which 
maximizes the transport, indicates the validity of the in- 
ternal mode mechanism also for the discrete case. 
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FIG. 13. Poincare sections for u> — 0.35, a — 0.1, 9 = 2, 
and N = 10 and different values of driving amplitudes and 
coupling constants. Panels (a)-(c) correspond to case k = 1. 
In panel (a), E l = E 2 = 0.182 (o), E x = E 2 = 0.18202 (+), 
Ei = E 2 = 0.18205 (dots). In panel (b), E ± = E 2 = 0.18207 
and in panel (c), E\ — E 2 — 0.1821. Panel (d) corresponds 
to k = 0.25 and Ei = E 2 = 0.097. 



IV. APPLICATIONS TO ARRAYS OF SMALL 
JOSEPHSON JUNCTIONS 

In this section, we discuss possible applications of 
above results to the case of an array of parallelly shunted 
and AC-biased small Josephson junctions (JJAs). This 
system is described by Eqs. (1) with u n correspond- 
ing to the phase difference of the wave functions of 



The voltage drop in an annular JJA (with one fluxon 
in it) is related to the average fluxon velocity (v) by the 
equation V = 2n(v)/N [35]. 

Figure 12 shows the IV characteristic of an annular 
JJA [dependence E\ (V)] for two different values of the 
coupling constant k. For k = 1 the characteristic is dom- 
inated by large phase locking steps associated to limit 
cycles with rotation numbers (k,l). In particular, the 
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phase locking resonances (k = I = 1) and (k = 2, 1 = 1) 
are very well pronounced. On the limit cycles of this type 
the voltage drop is equal to V = ku/(Nl) and for each 
lattice site n we have u n {t + IT) = u n {t) + 2nk (notice 
that for the phase locking dynamics, NV/uj corresponds 
to the ratio of the rotation numbers k/l). 

For smaller values of k (e.g., n — 0.25) directed fluxon 
motion becomes diffusive almost for all values of the driv- 
ing amplitudes. This behavior corresponds to the irregu- 
lar (non-vertical) parts of the IV curve seen for n — 0.25. 
Notice that in this case, only very few regular regimes are 
observed: pinning regions and very small phase locking 
steps with k = ±1, 1 = 1. An idea of complexity of the 
dynamics, resulting in this case, can be obtained from 
the Poincare sections of one of the junctions, as reported 
in Fig. 13 for different values of the driver amplitudes. 
The Poincare section is produced by plotting the values 
of the pair of the dynamical variables {u n {t),u n {t)} for 
n = 2 after the time intervals NT. The phase locking 
regime with the rotation numbers k = I = 1 is identified 
as the fixed point represented by the circle (o) in panel 
(a) of Fig. 13 and it corresponds to the fluxon which re- 
turns on the initial site (junction) after one trip around 
the array during the time NT. As the driving amplitude 
increases, the fixed point bifurcates into a period-two- 
orbit corresponding to the two crosses in the Poincare 
section shown in Fig. 13(a) (the kink dynamics repro- 
duces itself after the time 2NT, making two trips around 
the array, and so on). The curves shown in panel (a) 
refer to the case E\ = E 2 = 0.18205 and correspond to 
quasiperiodic motion. In panels (b) and (c) of Fig. 13 we 
report the Poincare sections of the period-doubling route 
to chaos, while panel (d) shows the case with n = 0.25 
for which chaos is fully developed and the fluxon motion 
is strongly diffusive. A similar behavior is found for the 
Poincare sections taken on other junctions in the array. 

Notice the presence of voltage sign reversals in the IV- 
characteristic in Fig. 12 generated by reversals of the 
direction of the motion of the kink. In particular, the 
sudden emergence of small intervals around E\ ~ 0.178 
and Ei ~ 0.11, associated with negative voltages at the 
resonances k = —1, I = 1 and k = —1, I = 3, respec- 
tively, are clearly visible. These phenomena, as discussed 
in the previous section, are due to attractor crises and 
have been observed also for other values n as well as 
for longer arrays (see also Fig. 8). Comparing the dy- 
namical behavior of arrays of different lengths (different 
number of junctions), we find that already for a small ar- 
ray of 10 junctions, the fluxon dynamics is qualitatively 
very similar to that obtained for the infinite chain. The 
quantitative difference, existing between the two cases, 
is probably due to the interaction of the fluxon with the 
radiated small amplitude waves (Joscphson plasmons), 
this being more pronounced in the discrete case due to 
smaller size of the system. 



V. CONCLUSIONS 

The ratchet phenomena induced by temporarily asym- 
metric zero mean fields on topological solitons (kinks 
and antikinks) of the discrete sine-Gordon equation have 
been investigated. In particular, we have studied the 
conditions for the occurrence of discrete soliton ratchets 
and the dependence of the average soliton velocity on 
the system parameters has been found. It has also be- 
ing shown that, in analogy with the continuous SG case, 
the unidirectional motion arises, when all symmetries of 
the system relating orbits with opposite velocities are 
broken. This condition can be achieved, using external 
biharmonic drivers of zero mean, which satisfy suitable 
inequalities. The main characteristic of discrete soliton 
ratchets, in comparison with the continuous case, is the 
existence of a depinning threshold (in the driver ampli- 
tude) related to the Peierls-Nabarro barrier. For driv- 
ing amplitudes smaller than the critical threshold, the 
kink performs periodic oscillations around its center of 
mass, while above this threshold, the kink directed mo- 
tion takes place. The ratchet dynamics are effectively 
achieved, when the kink motion is phase locked to the 
external driver. Besides phase locking regimes, we have 
shown that the transport is possible (but less effective) 
also in the presence of chaos and intcrmittency. From this 
point of view, the phenomenon appears to be much more 
complicated than the corresponding one in the continu- 
ous SG system. The dependence of the mean velocity 
of the kink on the amplitude, phase, and frequency of 
the biharmonic driver, as well as on the damping and 
on the coupling constant, have extensively been investi- 
gated. We have also found that the dependence of the 
velocity on the system parameters is much more compli- 
cated in the discrete model, resembling in most cases a 
devil's staircase. Some qualitative analogy with the con- 
tinuous case, however, remains. In particular, the valid- 
ity of the internal mode mechanism for discrete ratchets 
has been demonstrated. 

Finally, the possibility to observe experimentally dis- 
crete soliton ratchets in a one-dimensional AC-biased 
annular array of small Joscphson junctions has been 
discussed in detail. From our results we predict that 
the kink ratchet dynamics induced by a biharmonic AC 
driver of zero mean leads to the formation of voltage steps 
in the IV characteristic of the array, in the absence of any 
DC bias. 
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